Modeling the distribution of the endangered Jemez Mountains salamander (Plethodon neomexicanus) in relation to geology, topography, and climate

Abstract The Jemez Mountains salamander (Plethodon neomexicanus; hereafter JMS) is an endangered salamander restricted to the Jemez Mountains in north‐central New Mexico, United States. This strictly terrestrial and lungless species requires moist surface conditions for activities such as mating and foraging. Threats to its current habitat include fire suppression and ensuing severe fires, changes in forest composition, habitat fragmentation, and climate change. Forest composition changes resulting from reduced fire frequency and increased tree density suggest that its current aboveground habitat does not mirror its historically successful habitat regime. However, because of its limited habitat area and underground behavior, we hypothesized that geology and topography might play a significant role in the current distribution of the salamander. We modeled the distribution of the JMS using a machine learning algorithm to assess how geology, topography, and climate variables influence its distribution. The best habitat suitability model indicates that geology type and maximum winter temperature (November to March) were most important in predicting the distribution of the salamander (23.5% and 50.3% permutation importance, respectively). Minimum winter temperature was also an important variable (21.4%), suggesting this also plays a role in salamander habitat. Our habitat suitability map reveals low uncertainty in model predictions, and we found slight discrepancies between the designated critical habitat and the most suitable areas for the JMS. Because geological features are important to its distribution, we recommend that geological and topographical data are considered, both during survey design and in the description of localities of JMS records once detected.


| INTRODUC TI ON
Assessing the distributional extent of taxa is essential for species that are endemic, rare, and have limited dispersal capabilities, all of which increase their risk of extinction (Chunco et al., 2013;Smith & Green, 2005). This is because an accurate evaluation of habitat suitability informs management and conservation decisions by helping to identify potential areas to survey and protect, saving time and resources (Ancillotto et al., 2020;Crawford et al., 2020). Species Distribution Models (SDMs) help to evaluate habitat suitability and are key to estimating risk to species by understanding potential vulnerabilities and what landscape and climate variables contribute most to their distribution (Wang et al., 2020). These tools can also be useful for predicting species' responses to future climate and environmental change (Beest et al., 2021;Pang et al., 2021), which is especially important for climatically constrained species, such as those that live in mountainous areas (Ali et al., 2021;Rahbek et al., 2019;Xenarios et al., 2019).
The Jemez Mountains salamander (Plethodon neomexicanus, hereafter referred to as JMS) is endemic to the Jemez Mountains in north-central New Mexico and found primarily around the flanks and rim of the Valles-Toledo caldera complex in mixed-conifer forests (Degenhardt et al., 1996). It is a relatively rare and strictly terrestrial salamander and is the only Plethodon species in New Mexico.
Like other members of the family Plethodontidae, the JMS is highly fossorial and lungless, requiring moist conditions for surface activity, such as foraging or mating. It was listed as endangered under the Endangered Species Act in 2013 because of large and severe wildfires due to extensive drought (U.S. Fish and Wildlife Service, 2013).
The current federally delineated critical habitat of the JMS consists mostly of mixed-conifer forests. Principal threats to JMS habitat include historical fire exclusion and suppression, severe wildland fires, forest composition and structure conversions, postfire rehabilitation, forest and fire management, roads, trails, habitat fragmentation, recreation, and disease (U.S. Fish and Wildlife Service, 2013).
Jemez Mountains salamanders spend long periods of time underground, presumably in voids caused by the local geology, plant roots, or other processes; however, little is known about their underground habitat requirements. They move very little and have an estimated home range of 8 m 2 (Ramotnik, 1988). Because of their limited movements and underground behavior, we hypothesize that the geology and topography might play a significant role in their current distribution. Additionally, a recent study assessed the habitat suitability of members of the genus Plethodon in the Pacific Northwest and found that the distributions of these species depend strongly on precipitation (Nottingham & Pelletier, 2021) but the study's climate focus means it did not evaluate whether geologic and/or soil conditions contributed to species distribution.
Maxent is a type of machine learning SDM that uses a maximum entropy probability distribution to contrast occurrence data (i.e., presence records) with environmental data, such as climate, soil type, and geology, and estimates a probability distribution that has the maximum entropy (i.e., that is most spread out, or uniform) given certain constraints. The constraints are that the expected values of each feature, such as a climate variable, must equal the average value at known occurrence points, or most common value for categorical variables (Phillips et al., 2006). Maxent is one of the best algorithms to calculate the suitability of landscape for a species when presence/absence data are not available (Elith et al., 2006;. This is particularly the case for animals such as the JMS, because its absence is difficult to confirm due to its cryptic habits and low detectability. Models built with presence-only data do not incorporate information on the frequency of occurrence, and therefore cannot accurately predict probability of presence (Guisan & Thuiller, 2005;MacKenzie et al., 2002). However, such models can be used to estimate an index of the suitability of landscape for a species, including the relative importance of different variables (Elith et al., 2006).
Here, we provided an assessment of to what extent the geology, topography, and certain climate variables in the Jemez Mountains and around the Pajarito fault system influence the distribution of the JMS. We first compiled all available records of the occurrence of the salamander. We constructed Maxent models of potential habitat suitability at fine scales (5 m), by including features of geology, topography (i.e., elevation, slope, and topographical complexity from a LiDAR-derived digital elevation model [DEM]), and seasonal summaries of climate. We evaluated the relative importance of these climatic, geologic, and LiDAR-derived features in models of habitat suitability and determined which were meaningful and important parameters for identifying JMS habitat.

| Study area
The Jemez Mountains, situated in north-central New Mexico, includes the federally designated critical habitat of the JMS (Figure 1) and encompasses the Valles caldera, a resurgent volcanic caldera (Smith & Bailey, 1968). This volcanic system has been active since ca. 23 Ma (Gardner & Goff, 1984;WoldeGabriel et al., 2003), with two large ash-flow tuff-producing, caldera-forming eruptions at 1.6 and 1.25 Ma (age from Phillips et al., 2007), respectively. The two large ash-flow tuff-producing eruptions created the Bandelier Tuff, which blankets all of the surrounding flanks of Valles caldera. Subsequent smaller rhyolitic eruptions continued until 50-72 ka

T A X O N O M Y C L A S S I F I C A T I O N
Conservation ecology; Landscape ecology (Gardner et al., 1986;Wolff et al., 2011;Zimmerer et al., 2016). The Jemez Mountains and Valles caldera are separated from the Pajarito Plateau in the east by the Pajarito fault system, a dominantly downto-the-east complex normal fault that is potentially seismogenic and presently defines the active western margin of the Española Basin of the Rio Grande rift (Gardner et al., 1999;Lewis et al., 2002Lewis et al., , 2009). Elevation within the critical habitat ranges between 2200 and 3100 m (U.S. Fish and Wildlife Service, 2013). The eastern portion of the JMS habitat was recently subject to multiple large fires, some with high intensity and high severity burns and as large as 63,400 ha.
The climate within the JMS critical habitat is continental and semi-arid, with precipitation dominated by convective storms between July and September. Salamanders rely on this summer precipitation and resulting moisture for above-surface activity (Degenhardt et al., 1996). Winter precipitation is highly variable year to year due to Pacific Ocean teleconnections (Sheppard et al., 2002). Total annual precipitation within the critical habitat averages 643 mm (1981PRISM Climate Group, 2019). Annual average minimum temperature is −2.4 C, whereas the average maximum temperature is 12.3 C (1981-2010; PRISM Climate Group, 2019).

| Compilation of JMS presence records
We contacted several colleagues at different agencies and universities for data on surveys and any records of Jemez Mountains salamanders. In particular, we relied on field notes of the late Charles W. Painter (CWP; New Mexico Department of Game and F I G U R E 1 (a) Geologic map and area of analysis for modeling habitat suitability of the Jemez Mountains salamander in north-central New Mexico. The red points are localities where salamanders were detected in surveys or collected as specimens. These represent the thinned presence points (n = 189) that were used in the modeling process. This area includes the U.S. Fish and Wildlife Service federally designated critical habitat extent (outlined in black). (b) Number of salamander occurrences according to geology map units. The geology map unit colors in (b) match those in (a). Numbers in bars represent the percentage of salamander occurrences. Numbers in parentheses indicate the percentage of that geological map unit in our study area. Only those map units with greater than 10 salamander occurrences are presented. Descriptions of the map units in B can be found in the Appendix A.
Fish herpetologist) and other data deposited with the Museum of Southwestern Biology (MSB) at the University of New Mexico in association with JMS specimens. We queried, processed, and used information from Los Alamos National Laboratory (LANL; Hathcock et al., 2017), New Mexico Natural Heritage (NMNH), and the Global Biodiversity Information Facility (GBIF) data aggregator. We were able to compile 655 records from GBIF representing 49 distinct localities and reconciled these data with 690 records from NMNH, of which 650 were those documenting presence of the salamander. All presence records were georeferenced.
The records of JMS span from the years 1949 to 2017 and at each locality the number of salamanders documented ranges from 0 to 110 (mean = 4.86, median = 2) (110 = type locality, from where the species was described).

| Geological, topographical, and climate data
We used ten variables in our distribution models [geologic: (1) unit classification based on 1:24,000 scale geologic maps produced by New Mexico Bureau of Geology and Mineral Resources, New Mexico Institute of Mining and Technology Goff et al., 2002Goff et al., , 2012Kelley et al., 2004;Kempter et al., 2002;Timmer et al., 2006), (2) distance to the boundary of mapped geologic contacts (Goff et al., 2011) within the Valles caldera region; topographic: (3) high-resolution elevation (10 m), (4) slope, (5) topographic characterization (i.e., curvature; change in slope, first derivative) from a LiDAR-derived digital elevation model (DEM); climatic: (6) total precipitation in summer, (7) total precipitation in winter, (8) maximum temperature in winter, (9) minimum temperature in winter, (10) minimum temperature in summer]. were downscaled), we downscaled the data to match our intended 5 m fine scale for our analysis. To downscale the coarse scale PRISM data to 5 m, we followed the methodology of Lee et al. (2014). We believe that climate data, because of its continuous nature, is least likely to introduce bias at such scales.

| Modeling JMS distribution
We modeled the current extent of suitable habitats for JMS using Maxent implemented in the ENMeval package (Kass et al., 2021; version 2.0.0) within the R statistical framework (R Core Development Team, 2021). We intentionally set the area of analysis to include the federally designated critical habitat for the JMS (colored region of Figure 1a). We used this entire extent during the modeling and habitat suitability mapping process in order to find suitable habitats within, and outside, the federally designated critical habitat. The very northern part (<10 km 2 ) of federally designated critical habitat is not included in the analyses due to a lack of fine-scale geological data needed for modeling. Because of the gridded nature of analyses and models, we used all point coordinates and associated errors of JMS records as pixels representing presence of JMS (e.g., coordinates with a 10 m error would mean that 9 pixels are treated as if JMS was present there). In cases of multiple records, we only used one set of coordinates to minimize sampling bias, which is especially important for correlative modeling (Phillips et al., 2009). In addition, we spatially thinned the occurrences to 100 m using the thin func- We tested 4 types/combinations of feature classes: L, H, LQ, and LQH. We also tested four regularization multipliers; 0.5, 1, 2, and 5, which resulted in 16 total models. The higher the regularization multiplier, the higher the penalty for models with higher numbers of variables; thus, larger values encourage models with fewer covariates, lowering overfitting. All environmental variables were continuous variables, except for geological classification, which was set as a categorical variable.
We initiated Maxent to randomly sample 10,000 background points within our entire extent and trained the models using k-fold cross-validation using the jackknife partitioning method. Although We selected the best model using a combination of the area under the curve of the receiver-operating characteristic (AUC), omission rates, and AICc values. We use a custom stepwise process to select models to avoid overfitting (Gorris et al., 2021). First, we selected only half the models with the lowest omission rates using the 10% omission rate of the training localities metric (or.10p.avg).
Omission rates greater than the expected 10% usually means that the model overfits the data (Muscarella et al., 2014). From these models, we selected models that had the lowest difference between the training and testing AUC (auc.diff.avg) (Gorris et al., 2021). These models were chosen based on being lower than the median value of auc.diff.avg across the remaining models, which was 0.102. After the above steps, we then chose the model with the lowest AICc value as the top model. We assessed the included predictor variables in the top model using the built-in permutation importance and percent contribution. Together, these metrics allow for the identification of important predictors in the model (Cobos et al., 2019;Gorris et al., 2021).
After selecting the best model, we assessed the uncertainty in our predictions using the feature class and regularization multiplier of the best model. For each of 10 replicates, we bootstrapped our presence data (n = 189) using 80% of our presence data. We used the difference between the minimum and maximum habitat suitability (i.e., the range) among the 10 bootstrapped replicates to show areas of lower and higher uncertainty in our models (Gorris et al., 2021;Romero-Alvarez et al., 2020). For all models, we used 10,000 background points and the jackknife method of cross-validation.
We created a habitat suitability map for the entire extent using the mean values of the 10 bootstrapped models. We used the predict function to create a habitat suitability map (occurrence intensities) for each of the bootstrapped models. We used the complementary log-log (cloglog) transformation to give us probabilities of presence. We used these transformed model outputs because they fall between 0 and 1. However, the term "probability of presence" is subject to a few assumptions about the sampling F I G U R E 2 (a) Habitat suitability map for JMS from the Maxent models. Colors indicate the mean habitat suitability from the 10 bootstrapped Maxent models using the parameters of the top model after model selection. The northern part of the range is not included in the map due to lack of fine-scale geological data available for modeling. Yellow colors indicate areas with high habitat suitability, while darker blue colors indicate areas with lower habitat suitability. The red outline is the federally designated critical habitat designated by the U.S. Fish and Wildlife Service. (b) Map of the study area depicting the uncertainty in habitat suitability for the Jemez Mountains Salamander. Colors indicate the range in maximum and minimum values in habitat suitability from the 10 bootstrapped Maxent models using the parameters of the top model after model selection.
scheme (Phillips et al., 2017). Therefore, we describe the outputs as "relative habitat suitability" with 0 being low habitat suitability, 0.5 being medium suitability, and 1.0 being high suitability (Gorris et al., 2021). The permutation importance and percent variable contribution are reported as the means from all 10 bootstrapped models. these coarse comparisons only suggest patterns with geology and TA B L E 1 Summary of the four top Maxent models that passed the omission rate and difference between training and test AUC thresholds (see Section 2). Included here are types of feature classes, regularization multipliers, AUC for training data, AICc values, the deviation from the best model (ΔAICc), and number of model parameters. The top model has a feature class of LQ, a regularization multiplier of 2, and 35 parameters.

| Modeling JMS distribution
We used Maxent to find the most important geological, topographical, and climate variables important for the distribution of the JMS.
The top model for JMS habitat suitability had a linear and quadratic (LQ) feature class, a regularization multiplier of 2, and 35 model parameters ( Table 1). The selection of this model was based on a custom set of thresholds in order to not overfit the data. The variable with the highest percent contribution was geological classification (45.6%) followed by maximum temperature during the winter months (26.3%; Table 2). The variable with the highest percent permutation importance is maximum temperature during the winter months (50.3%) followed by geological classification (23.5%;

| DISCUSS ION
Most Plethodontid salamanders have limited ranges and many species need protection due to habitat vulnerability (Milanovich et al., 2010). The Jemez Mountains salamander is endemic to New Mexico and more specifically to the flanks of the Valles caldera in mixed-conifer forests (Degenhardt et al., 1996). Further threats to this federally listed endangered species include declining or changing forest cover, changing fire regimes resulting in less frequent but more severe fire, increases in temperatures of soil, and associated evaporation, and changes in precipitation patterns (U.S. Fish and Wildlife Service, 2013 In some portions of the study area, topographic complexity and geology can be correlated, suggesting a potential importance of such areas for the JMS. For example, a geologically related element, but one that is not incorporated as a model input, There were several geologic map units in which salamander occurrences were greatest. Qbt (Bandelier Tuff, Tshirege Member) had the most salamander occurrences (29.6%). This is consistent with this geology type being the most frequent in our study area (28.81%; Figure 2b). However, the geology types Qbo, Tpa, Qcbt, and Ttcg also had high salamander occurrences. For these areas, salamanders are more numerous than the frequency of these geology types in our study area (ratio of salamanders to geology type >3), suggesting that salamanders congregate in these geology types since they are disproportionately inhabited by salamanders. The occurrence of JMS in certain geology types may represent a correlation with certain soil characteristics. Geology is often a coarse surrogate for soil characteristics, which are often an important component of the habitat of certain salamander species. In addition to the geological and topographical variables considered here, future work should consider additional variables, such as more accurate soil composition, pH, and moisture retention (Nottingham & Pelletier, 2021 Hathcock: Conceptualization (equal); investigation (equal); writingreview and editing (equal).

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.  (Spell & Harrison, 1993); magnetic polarity normal; maximum exposed thickness at least 510 m

Qvsm2
Quaternary South Mountain Member-Flow-banded, massive to slightly vesicular porphyritic rhyolite lavas containing abundant phenocrysts of sanidine, plagioclase, quartz, biotite, hornblende, and clinopyroxene in a pale-gray, perlitic to white, devitrified groundmass; apparently consists of four flow units based on morphology (youngest to oldest Qvsm4 to Qvsm1); fills paleocanyon in southern moat of Valles caldera; 40Ar/39Ar age of Qvsm2 is 0.52 ± 0.01 Ma (Spell & Harrison, 1993); maximum exposed thickness is at least 450 m. Includes Cerro La Jara rhyolite (Qvlj), a small dome of flow-banded, massive to slightly vesicular porphyritic lava; 40Ar/39Ar age is 0.53 ± 0.01 Ma (Spell & Harrison, 1993); magnetic polarity normal; maximum exposed thickness about 75 m Tbh Tertiary Bearhead Rhyolite-Dikes, plugs, and flows of aphyric to slightly porphyritic, devitrified to completely silicified rhyolite containing sparse phenocrysts of quartz, sanidine, plagioclase, biotite, opaque oxides ± hornblende; locally shows pervasive hydrothermal alteration consisting of quartz, chalcedony and/or opal, illite, Fe-and Mn-oxides, pyrite, and possibly other sulfides, alunite, jarosite, and gypsum; 40Ar/39Ar ages on widely separated samples range from 4.81 to 7.83 Ma (Justet & Spell, 2001;Kempter et al., 2007) maximum observed thickness about 100 m Tpa Tertiary Two-pyroxene andesite, undivided-Domes, flows, flow breccia, spatter deposits, and scoria of andesite from multiple sources; vents are widely scattered; individual units are slightly porphyritic to very porphyritic containing phenocrysts of plagioclase, orthopyroxene, and clinopyroxene; alteration varies from slight to intense consisting of silica, calcite, Fe-oxides, clay ± chlorite ± zeolite ± pyrite ± epidote; 40Ar/39Ar ages in western and southern map area range from 8.2 to 9.4 Ma (Justet, 2003); maximum exposed thickness about 150 m Tpb Tertiary Olivine basalt and basaltic andesite, undivided-Flows, flow breccia, spatter deposits, and scoria of basalt and subordinate basaltic andesite from multiple vents; most units are slightly porphyritic containing phenocrysts of olivine, plagioclase ± clinopyroxene; displays variable amounts of hydrothermal alteration consisting of silica, calcite, Fe-oxides, clay ± zeolite ± chlorite ± epidote ± pyrite; 40Ar/39Ar ages from western and southern map areas range from 8.88 to 9.45 Ma (Justet, 2003); maximum exposed thickness about 150 m Tpbhd Tertiary Porphyritic biotite, hornblende dacite-Extensive dome and flow complex filling paleocanyon south of Rabbit Mountain; contains large phenocrysts of plagioclase, plus biotite, hornblende, orthopyroxene, and clinopyroxene; contains vesiculated enclaves of plagioclase, pyroxene ± hornblende ± biotite as large as 30 cm in diameter; hydrothermally altered to clay, silica, calcite, Fe-oxides, chlorite ± epidote; 40Ar/39Ar age is 8.66 ± 0.22 Ma; exposed thickness at least 275 m Tpv Tertiary Volcaniclastic member (Pliocene? to Miocene)-Tpv is conglomeratic sandstone and sandy conglomerate locally containing cinder deposits, pyroclastic-fall deposits, and lava flows too small or thin to map; unit has accumulated in small basins, topographic lows, and paleocanyons; contemporaneous with eruption of lavas of the Paliza Canyon Formation; upper part of unit may be correlative with oldest Cochiti Formation (Smith & Lavine, 1996); maximum exposed thickness about 70 m Tpvs Tertiary volcaniclastic sandstone-Tpvs is a moderately to well-sorted volcaniclastic sandstone that is brick red to tan and contains mostly volcanic fragments, feldspar, mafic minerals, and minor quartz; present between lava flow contacts in isolated locations throughout southeastern part of the map area; mapped only where laterally extensive and at least 3 m thick Ttcg Tertiary Cerro Grande dacite-Extensive dome and flow complex of massive to sheeted, porphyritic lava containing phenocrysts of plagioclase, hypersthene, and (typically) conspicuous hornblende; the latter two phases commonly show oxidized rims that may be difficult to see in hand sample; ages on widely separated samples range from 2.88 to 3.35 Ma (Broxton et al., 2007); maximum exposed thickness is about 750 m Ttpm Tertiary Pajarito Mountain dacite-Dome and flow complex of massive to sheeted, porphyritic lava containing phenocrysts of plagioclase, hypersthene, clinopyroxene, and opaque oxides in a devitrified groundmass; 40Ar/39Ar ages on geographically separated samples range from 2.93 to 3.09 Ma (Broxton et al., 2007); maximum exposed thickness is about 365 m Ttrc Tertiary Rendija Canyon rhyodacite-Dome and flow complex of massive to sheeted, highly porphyritic lavas with phenocrysts of quartz, sanidine, plagioclase, hornblende, and biotite; 40Ar/39Ar ages on widely separated samples range from 3.50 to 5.36 Ma (Broxton et al., 2007); maximum exposed thickness approximately 500 m